Particle simulation study of electron heating by 
counterstreaming ion beams ahead of supernova 
remnant shocks 

o 

<N 

Ph ■ ME Dieckmann 1 , A Bret 2,3 , G Sarri 1 , E Perez Alvaro 3 , I 

Kourakis 1 and M Borghesi 1 

1. Queen's University Belfast, Ctr Plasma Phys, Belfast BT7 INN, UK 

2. Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA 

3. ETSI Industrials, Universidad de Castilla-La Mancha, 13071 Ciudad Real, Spain 
and Instituto de Investigaciones Energeticas y Aplicaciones Industriales, Campus 

| Universitario de Ciudad Real, 13071 Ciudad Real, Spain 



^p H ' Abstract. The growth and saturation of Buneman-type instabilities is examined 

with a particle-in-cell (PIC) simulation for parameters that are representative for 
| the foreshock region of fast supernova remnant (SNR) shocks. A dense ion beam 

' and the electrons correspond to the upstream plasma and a fast ion beam to the 

shock-reflected ions. The purpose of the 2D simulation is to identify the nonlinear 
saturation mechanisms, the electron heating and potential secondary instabilities that 
^ ■ arise from anisotropic electron heating and result in the growth of magnetic fields. We 

confirm that the instabilities between both ion beams and the electrons saturate by the 
' formation of phase space holes by the beam-aligned modes. The slower oblique modes 

, accelerate some electrons, but they can not heat up the electrons significantly before 

\Q | they are trapped by the faster beam-aligned modes. Two circular electron velocity 

■ distributions develop, which are centred around the velocity of each ion beam. They 

develop due to the scattering of the electrons by the electrostatic wave potentials. The 
growth of magnetic fields is observed, but their amplitude remains low. 



J3 ; PACS numbers: 52.35. Qz, 52.50.Gj, 52.65.Rr 
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1. Introduction 

The thermalisation of shock-reflected ion beams plays an important role in solar system 
and astrophysical collision-less plasma. Energetic beams, which consist of the ions that 
were reflected by a plasma shock, outrun the shock and interact with the upstream 
plasma. This upstream plasma is the interstellar medium (ISM) for supernova remnant 
(SNR) shocks and the solar wind for most solar system shocks. The ion beam is a source 
of free energy, which is released through wave instabilities. Such instabilities have been 
observed ahead of solar system shocks pp. They are also important mechanisms for 
particle acceleration at supernova remnant (SNR) shocks [2-9]. 

The Buneman instability (TOj [11] develops between one ion beam and one electron 
beam of equal density, which are both spatially uniform, unmagnetized and drift relative 
to each other. Such a plasma is not current neutral. A plasma with no net charge and 
no net current can be composed of two counter-streaming ion beams and one electron 
species, if the total charge density of the ions is that of the electrons and if the partial 
currents of the three beams cancel out each other. This can be an appropriate description 
of the plasma upstream of a shock. One ion beam is composed of the shock-reflected 
ions and the second ion beam and the electrons are provided by the upstream plasma, 
into which the shock expands. Such distributions are observed ahead of the Earth's bow 
shock pQ. We refer with Buneman- type instability (BTI) to an electrostatic instability, 
which involves an ion beam with a density below that of the electrons. 

A BTI develops, if the drift speed between the ion beam and the electrons exceeds 
the electron thermal speed. Thermal damping effects due to the ions can be neglected, 
as long as they are reasonably cold. If the shock-reflected ion beam is much faster than 
the electron thermal speed and sufficiently dense, then its current has to be cancelled out 
by a high drift speed between the background electrons and ions. This drift speed may 
exceed the electron thermal speed and result in a second BTI. The speed of SNR shocks 
is a few per cent of the light speed c, which exceeds by far the electron thermal speed in 
the ISM. We may expect that the shock-reflected ion beam and the counterstreaming 
ion beam are both sufficiently fast in the electron rest frame to render the plasma 
unstable. Such plasmas have been examined widely in the past with particle-in-cell 
(PIC) simulations [12j. Previous studies have addressed unmagnetized [131 El 03] 
and magnetized plasmas [To] ITT] with one-dimensional PIC simulations, which can not 
capture the multi-dimensional nature of the wave fields. More recently, the BTI has 
been examined with two-dimensional PIC simulations in unmagnetized and magnetized 
plasma. The electrostatic simulation in Ref. [18] has modeled the interaction of one fast 
ion beam with cool electrons. References [HI EH] investigated ion beams that stream 
with nonrelativistic and mildly relativistic speeds across an orthogonal magnetic field. 

We examine here with an electromagnetic 2D PIC simulation how a plasma 
thermalises, which consists of two counterstreaming ion beams with unequal densities 
that move through an unmagnetized electron plasma. The drift speeds between the 
electrons and each of the ion beams exceed the initial electron thermal speed by factors 
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of 3.4 and 17 and two BTI's develop. The faster beam has a speed that is comparable 
to the ones modeled with PIC simulations in Ref. [18]. The waves driven by both BTI's 
interact nonlinearly with the electrons, once their amplitude is adequate. Trapping by an 
electrostatic wave accelerates electrons along the wave vector and a thermal anisotropy 
in the electron's velocity distribution develops as we show here. It is the purpose of 
our study to assess, if a thermal anisotropy-driven Weibel instability (TAWI) [21-29] 
is triggered by anisotropic electron heating. Such a secondary instability can not be 
resolved by the electrostatic PIC code used in Ref. [T8] . 

The initial conditions, which we consider here, are representative for a volume 
element in the foreshock, that is so small that the flow speed of the upstream plasma 
and of the shock-reflected ion beam are constant inside. However, ion beam-driven 
instabilities can also be triggered by the return current [30] of the cosmic rays. Similar 
instabilities may thus occur in a much larger volume ahead of the shock. The TAWI 
could provide [31] the seed magnetic field, which speeds up the cosmic ray-driven 
instabilities that magnetize SNR shocks j32j [33]. The inital magnetic amplitude is set 
to zero in our simulation, which simplifies the detection and interpretation of magnetic 
fields that are driven by the BTI through the subsequent TAWI. This choice can be 
justified by the weak magnetic fields ahead of SNR shocks, which yield an electron 
cyclotron frequency well below the BTI's frequency in the electron rest frame. At the 
same time, the magnetic field close to SNR shocks is thought to be strong enough to 
reflect enough ions to form a dense beam, provided that the shock is perpendicular. 

Our results are as follows. The interaction of the background electrons with the slow 
and dense ion beam from the background plasma yields the faster-growing instability, 
as we expect from the solution of the linear dispersion relation. The electrostatic waves 
grow and saturate through the formation of electron phase space tubes [34-39]. The 
spectrum of the unstable waves is not unidirectional [18] and the phase space tubes are 
spatially bounded. Their extent orthogonal to their wave vector is a few times their 
wave length. Trapping accelerates the electrons into the direction of the propagating 
waves. The electron thermal energy along the beam direction reaches twice the value 
of the electron thermal energy orthogonal to the beam direction and this anisotropy 
results in the growth of magnetic fields. The magnetic fields do, however, only reach an 
energy density that is 2-3 orders of magnitude below that of the electric fields and their 
strength does not exceed that of the interstellar magnetic fields. 

The electrons are initially accelerated by trapping to a speed, which is about twice 
as high as the background ion speed in the electron rest frame. After this initial 
acceleration, the electrons are scattered by the waves. The reflections are elastic, because 
the electrostatic wave potential is almost time-stationary in the reference frame of the 
bulk ions, and they change the momentum direction of the electrons. Repeated collisions 
between the electrons and the phase space tubes redistribute the electrons into a circular 
interval in velocity space, which is centered at the phase speed of the waves. This 
mechanism differs subtly from that proposed in Ref. [18], which has attributed the 
crescent formation to the low phase speed of the oblique modes. 
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The plasma is no longer free of current after the saturation of the first BTI. The 
electrons have been decelerated in the rest frame of the background ions by their 
interaction with the waves while the ions with their high inertia have only been weakly 
affected. A macroscopic beam-aligned electrostatic field grows in response to this current 
on scales much larger than the wave length of the BTI [40J. This field accelerates the 
electrons nonresonantly such that the mean electron speed remains close to the one prior 
to the wave saturation, thereby reducing the net current. 

Eventually the electrons start to interact with the faster ion beam. Electrons are 
accelerated primarily by trapping, but the Landau damping of the oblique modes results 
in a separate population of hot electrons. A crescent distribution develops, but it is less 
pronounced than the one found in Ref. [18] . The electron heating by the first BTI 
and the electron acceleration by the macroscopic electric field have changed the electron 
distribution function to a non-Maxwellian one and they have decreased the ratio between 
beam speed and electron thermal speed to a value well below that used in Ref. [18J, 
which may explain the different plasma evolution. The hotter ion beam we consider here 
will also modify the wave spectrum of the BTI and the nonlinear response of the beam 
ions to the electric fields, which may also contribute to the different interaction of the 
electrons with the waves driven by the second BTI. An energetic circular distribution 
in velocity space gradually develops after the second BTI has saturated by electron 
trapping. The cause is again the electron scattering by the electrostatic structures that 
move now with the mean velocity of the fast ion beam. No TAWI develops in this 
case during the simulation time, even though the electrons continue to show a thermal 
anisotropy. The absence of any significant magnetic field growth during the growth and 
saturation of the two BTIs suggests that the thermal anisotropy created by the nonlinear 
interaction between electrostatic waves and electrons is not capable of magnetising the 
foreshock of supernova remnant shocks, at least not for our initial conditions. 

The structure of the paper is as follows. Section 2 discusses the initial conditions 
and the PIC simulation code. Section 3 presents the simulation results. Section 4 is the 
discussion. 

2. The initial conditions and the particle-in-cell simulation code 

2.1. Initial conditions 

We consider here three spatially uniform plasma species with non-relativistic Maxwellian 
velocity distributions. The electrons have the density no, the plasma frequency 
w p — { n o£ 2 /^e^o) 1 ^ 2 and the temperature 10 eV or the thermal speed v te = 1.325 x 10 6 
m/s. They move at the speed v e = 4.5 x 10 6 m/s along increasing values of x, which 
gives v e = 3Av te . We consider for computational reasons a reduced mass ratio between 
the ions and electrons of mi/m e = 10 3 . The bulk ions are at rest in the simulation 
box. Their density and temperature are = 5uq/6 and 5.5 eV. The bulk ions would 
correspond to the relatively cool upstream plasma ahead of a SNR shock. The ion beam 
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Figure 1. (Color online) The beam distribution is displayed in panel (a). Panel 
(b) shows the numerical solution of the linear dispersion relation as a function of the 
beam aligned wave number k x c/uj p and the perpendicular wave number k y c/aj p . The 
color corresponds to the exponential growth rate of the wave, expressed in units of 
the electron plasma frequency uj v . The unstable wave branch at large k x arises from 
the interaction between background ions and electrons and that at low k x from the 
interaction between electrons and beam ions. 

has the density n& = no/6 and temperature 10 keV. Its mean speed is Vb = 6v e . The 
ion beam is much hotter than the bulk ions. Plasma shocks are neither perfectly planar 
nor elastic reflectors, so the momentum change by the reflection differs for each ion [7j, 
which causes a rise of the temperature. The velocity distribution of the three species is 
displayed in Fig. Ufa). The plasma is charge neutral by + = n and it is initially 
free of current by n v e = n^Vb. 

Both BTI's result in electrostatic waves that are practically stationary in the rest 
frame of the respective ion beam. The real parts of the wave frequencies u u i,co U 2 and 
the corresponding wave numbers k ul , k u2 along x can be estimated in the cold plasma 
approximation with the dispersion relation 1 — u 2 /u 2 — u 2 / (co — k x vo) 2 = valid for one 
ion beam with plasma frequency Ub and the electrons with the drift speed vq. The bulk 
ions and the electrons drive waves with the frequency uo u \ <C 1 and the wave number 
k u ic/up ~ c/v e or k u ic/up ~ 67. Since the reference frame of the bulk ions is that of 
the simulation box, this wave will be practically stationary in the latter. 

The instability between the fast ion beam and the electron beam drives waves 
with k u2 c/u p ~ c/(vb — v e ) or k u2 c/u p ~ 13. Their Doppler shifted frequency is 
u U 2 ~ k u2 Vb or u u2 ~ 1.2 in the box frame of reference. The growth rate of a BTI 

1/3 

is LOi ~ (3 \/3lo 2 uj p / '16) [13J and the growth rate ratio between both BTI's is thus 
^ii/^i2 = (ni/nb) 1 ^ 3 . The exponential growth rate of the BTI driven by the bulk ions is 
1.8 times that of the BTI driven by the fast ion beam. 

In the cold plasma limit, this growth rate is independent of the wavenumber k y 
orthogonal to the beam velocity vector. A more accurate growth rate map is obtained 
from the solution of the linear dispersion relation that takes into account thermal effects 
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and is valid for waves in the simulation plane with k 2 = k 2 + k 2 

(u 2 e xx - k 2 y c 2 ^ {uj 2 e yy - k 2 x c 2 ^ - (u 2 e yx + k x k y c 2 ^ = 0, (1) 
and with the elements of the dielectric tensor 

o 2 k-(^] 

e a ,(k, W ) = ^ + E^f/ v a d Jld\ + E 4 / W P \ V d\ (2) 
j u J ovp j to 2 J ijj — k • v 

where u P j and /? are the plasma frequency and equilibrium distribution function of 
the j th species, respectively (See section 2 in |JT]). We approximate f® by waterbag 
distributions with the same densities and drift velocities as the plasma species we 
introduce in our simulation. The thermal width of each waterbag distribution equals the 
thermal speed of the corresponding Maxwellian distribution. The waterbag distribution 
is equivalent to a warm fluid model, and a good approximation to a Maxwellian, if the 
beam speed is large compared to the thermal spreads |42j . The solution of the linear 
dispersion relation is shown in Fig. QJb), which displays the imaginary part of the wave 
frequency (growth rate) as a function of the beam-aligned wave number k x and of the 
perpendicular wave number k y . The colour corresponds to the growth rate, expressed 
in units of the electron plasma frequency u p . The most unstable wave numbers in this 
growth rate map match the ones from the cold plasma approximation for k y ~ 0. Larger 
values of k y result in increasing values of k x and in lower growth rates. 

2.2. The simulation code and the diagnostics 

Our 2D3V electromagnetic and relativistic PIC simulation code (13] resolves the x, y 
plane and the particle positions in this plane. It updates all three components of the 
magnetic field B, of the electric field E, of the current J and of the relativistic particle 
momentum p. The code preserves V • E = p/eo and V • B = to round-off precision 
and evolves the electromagnetic fields in time through 

<9B <9E 
V x E = -— , V x B = mo-^- + /A)J- (3) 

An ensemble of computational particles (CPs) is followed in time. Each CP with 
the index i of the species j in the simulation has the mass and charge irtj and qj, 
the momentum = m^TiVj and the position x« = (xi,yi). Each CP corresponds to 
a volume element of the phase space distribution of the particle species it represents 
and the ratio qj/mj must be equal to that of the corresponding physical particle. The 
position and momentum of each CP are updated according to 

^ i = v 4 ,^ = g J (E[x,]+v l xB[x,]). (4) 

The simulation box has the size L x x L y = 3c/u p x Ac/u p , which is resolved by a grid 
with N x x N y = 580 x 760 quadratic cells with the sidelength A x w 1.2 Ad, where 
Ad is the electron Debye length. The boundary conditions are periodic, which limits 
the wavenumber spectrum to k x c/u p = 27ra/3 and k y c/u p = 2ixb/A with the integers 
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1 < a < N x /2 and 1 < b < N y /2. The short unstable wave has k u \L x « 200 and about 
32 wave lengths are resolved. The long unstable wave has k U 2L x « 39 and about 6 wave 
periods are resolved. Each of the plasma species is represented by 320 CPs per cell. The 
simulation time T M u p = 1200 is resolved by the time steps A t u p = 2.5 x 10~ 3 . 

We compute the following quantities to monitor the plasma evolution. The electric 
field energies D ex (t) = (2e ) 1 J E x (x,y,t)dx dy and D ey (t) = (2e ) _1 / E y (x,y,t)dx dy 
show the time-evolution of the BTI's. The energy D hz {t) = (2/i )~ J B 2 (x,y,t)dx dy 
reveals if an electron thermal anisotropy A ^ with A = (v x )/(Vy) — 1, (v x ) = 
J2 S ( v x,s ~ v e) an d (Vy) = J2 S v y, s resu lt s m a magnetic field growth. The summation 
is here over all computational electrons. The electron mean speed along x and y varies 
in the simulation only by about 10~ 2 v te and can thus be taken to be constant. The 
kinetic energies of the electrons and of the bulk ions are computed in our relativistic 
code as K e {t) = m^c 2 X) s (r s — 1) and Ki(t) = rriic 2 ^2 S (T S — 1), where we and mi are 
the non-relativistic masses of the computational electrons and ions. The sum is over all 
CPs of the respective species and T s the Lorentz factor of the s th CP. All energies are 
computed in the rest frame of the bulk ions and they are normalized by K e (0) = E^o- 
The mean electric field in the simulation box is (E x ) = (A^Ay) -1 J2iJi J2f=i E x (i, j, t), 
where E x (i,j, t) is the electric field at the grid cell with the indices i and j at the time t. 
A low pass filter removes the oscillations with a frequency exceeding Uf = 0.8u p , which 
reveals the trend of the curve. 

In what follows, the velocities v and positions x are normalized to v e and to the 
electron skin depth X s = c/u p . The time is multiplied with u p . Frequencies are given in 
units of Up and the wave numbers are given as kc/u p . The electric and magnetic fields 
are expressed as e^/ui p cm e and eH/cu p m e . 

3. Numerical simulation results 

Electrostatic waves are polarized parallel to their wave vector k. Their spectrum 
contains oblique modes and E x and E y will both grow. Figure |2] shows the time evolution 
of their energies. An exponential growth of the electric field energy is observed in the 
time interval 50 < t < 150 with D ex w 10D ey . Fitting exp (2cUji) to the energy densities 
D ex and D ey gives us the amplitude's growth rate Ui ~ 0.035. This growth rate is well 
below the peak un ~ 0.058 in Fig. [T](b). D ex and D ey are integrated over the entire 
simulation box, which averages the field energies over all unstable modes and not just 
over the fastest-growing ones. The saturation of D ex and D ey coincides with the growth 
of K e and Ki, which evidences the saturation of the BTI driven by these two species. 

Figure [2(b) reveals that d t {E x ) < during 100 < t < 150, which implies by 
Ampere's law a box-averaged {k x = 0) net current j x > and thus an electron 
deceleration in the reference frame of the simulation box. The net current is tied to 
a change of the mean speed of the electrons in response to the BTI between the bulk 
ions and electrons and their re-distribution in phase space by its saturation. This current 
is not compensated by a change in the ion mean speed due to the large ion inertia. The 
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Figure 2. (Colour online) The electric D ex (upper curve), D ey and the magnetic 
Db z (lower curve) energies are displayed in panel (a). The dashed blue lines are 
exponential fits with the same growth rate. Panel (b) shows the box-averaged electric 
field amplitude (E x ). The electron energy K e (t) is shown in panel (c) and that of the 
bulk ions in panel (d). All energies are normalized to Eko- 



electrons are accelerated by the (E x ) < 0, so that the electron mean speed along x 
remains close to v e . The required acceleration of electrons by (E x ) is small. Even the 
peak electric field (E x ) ~ —4 x 10 -5 can change the electron speed by only A„ f» 3 x 10 -3 
during a time interval A T 1. A rapid switch to d t {E x ) > takes place at t m 150, 
which coincides with the onset of the growth of the electron and bulk ion energies and 
is thus related to the nonlinear saturation of the BTI between these species. 

All displayed field components and the electron energy continue to grow after 
t ~ 150. D ex and D ey grow by another order of magnitude until t 500 when they 
reach their maximum. We show below that this growth is caused by the BTI between 
the fast ion beam and the electrons. The growth rates of D ex and D ey are well below 
that of the initial growth phase prior to t ~ 150 and D ex is not proportional to D ey . We 
observe a d t (E x ) > and thus a change of the sign of the net current to positive values 
after t ps 350. Electrons start to interact with the waves driven by the second BTI and 
are accelerated to positive v x . An exponential growth of D bz between 200 < t < 400 is 
observed and it remains approximately constant after this time. The magnetic energy 
remains 2-3 orders of magnitude below that of the electric field components, which 
implies an essentially electrostatic plasma dynamics. The electron energy K e , which 
takes into account their bulk flow, grows practically linearly up to a value QEko i n the 
interval between 150 < t < 1000. 

3.1. Saturation of the instability between electrons and bulk ions 

Figure [3] shows the electric fields in a part of the simulation box at t = 150. The 
structures in E x and E y belong to waves with wave vectors that are in some cases tilted 
with respect to the beam velocity vector. Their amplitude is two orders of magnitude 
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(a) X-position (b) X-position 



Figure 3. (Color online) The electric field distributions in a subinterval of the 
simulation box at t = 150: E x is shown in panel (a) and E y in panel (b). 




(b) X-position (c) \ 



Figure 4. (Color online) Electron phase space distributions at t = 150: Panel (a) and 

(b) show phase space projections onto the x, v x plane. The distribution is integrated 
over < y < 0.05 in (a) and over 1.55 < y < 1.6 in (b). The electron distribution 
has been integrated over all x, y,v z for the projection onto the v x ,v y plane in panel 

(c) . The vertical line corresponds to v x = v e . All distributions are normalized to their 
respective maxima and the colour scale is 10-logarithmic. 

larger than the mean electric field (E x ) in Fig. E(b). Such wave fields have also been 
observed previously in electrostatic PIC simulations [18]. The E x and E y are projections 
of the obliquely oriented electric field onto the x and y axes, explaining why D ex oc D ey 
until t = 150. The wavelength of the structures is A ~ 0.06 or k = (k^ + ky) 1 ^ 2 ~ 100. 
This value is in line with the wave numbers of the oblique modes in Fig. [U^b). 

Various projections of the resolved electron phase space distribution f(x, y, v x , v y , v z ) 
at t = 150 are shown in Fig. HJ The phase space projection on the x, v x plane integrates 
this distribution over all v y , v z and over ten grid cells along y. Figures H^a) and (b) inte- 
grate the distribution over < y < 0.05 and over 1.55 < y < 1.6, respectively. Figures 
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Figure 5. (Color online) Evolution of the electron velocity distribution: Panel (a) 
shows the electron distribution at t = 200 and panel (b) shows it at t = 300. Both 
distributions are normalized to the peak value at t = 150 and the colour scale is 
10-logarithmic. The black curves are the contour lines -1.24, -1.23, -1.22, -1.21, -1.20. 



S(a) and (b) reveal electron phase space holes, which arise when electrostatic waves sat- 
urate [31]. The spatial width of each phase space hole corresponds to one wave period 
of the electrostatic waves driven by the BTI. The exponential growth phase of D ex and 
D ey thus ends at t — 150 as usual due to the trapping of electrons. The electrons are 
trapped by a potential, which is practically stationary in the rest frame of the bulk ions. 
This is evidenced by the supplementary movie 1, which shows a time animation of Fig. 
@|a). The electrons gyrate around v x = at early times. The mean speed of the trapped 
electrons inside of a phase space hole is less than v e , which can be seen in Fig. HI A 
significant fraction of the electrons has been slowed down at t = 150, which reduces the 
modulus of the electron current. The electron current is negative because v e > 0. The 
ions hardly react to the electric field because rrij ^> m e and they can not compensate the 
change of the electron current. We obtain a macroscopic box-averaged (j x ) > 0, which 
explains why d t (E x ) < in Fig. MJd). The mean electric field (E x ) < in Fig. [2(b) 
accelerates all electrons in the positive direction, which will reduce (j x ). The maximum 
of the phase space density in Fig. @](c) is thus found at v x > 1. The electrons in this 
core population have not been trapped as we can see from Figs HIa,b) and from movie 
1, but they have all gained speed along x. The trapped electrons in Fig. @Jc) are found 
at v x < 1. The electron distribution is stretched out along v y , because the electrons are 
trapped by oblique waves (See Fig. [3]). 

The electrons are accelerated by (E x ) < (Fig. E^b)) into the positive x-direction 
and they are heated by their interaction with the localized electric field structures (Fig. 
EJ). Figure [5] demonstrates this by a comparison of the electron velocity distributions at 
the times t = 200 and t = 300, when the modulus of (E x ) is largest. The velocity that 
corresponds to the maximum of the electron phase space density is found at v x ~ 1.2 
at t = 200 and at v x ~ 1.6 at t — 300. The distribution widens in velocity space, which 
corresponds to an increasing thermal speed, and it adopts an increasingly circular shape. 
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(a) X-position (b) X-position (c) Time 



Figure 6. (Color online) The B z component of the magnetic field at the times t = 150 
(a) and t — 400 (b). The full simulation grid is displayed. The thermal anisotropy 
A is shown as a function of time w p t in (c) and the times t — 150 and t = 400 are 
overplotted as vertical dashed lines. 



The centre of the hot circular distribution is at v x ~ 0.3 at t — 300. We can not observe 
a density distribution in Fig. [5] that resembles the crescent found in Ref. |18j . 

The electron distribution in Fig. |4]^c) can be subdivided into a thermal core 
population centred at v x ~ 1.2 and v y with a high density and a thermally 
anisotropic hot electron population with a lower density. The low-density electron 
population is elongated along v x in Fig. |5]^a) and it has become almost circular at 
t = 300 in Fig. |5^b). It is thus fair to assume that this thermally anisotropic electron 
population is responsible for the growth of Df, z in Fig. |2|a) and that the magnetic field 
growth ceases once the electron distribution is isotropic in velocity space. A thermally 
anisotropic electron distribution drives waves with a wave vector, which is aligned with 
the direction along which the electrons are cool [2T]. We thus expect a wave vector of 
the magnetowaves that is almost parallel to the y-direction. The magnetic instability is 
driven by currents in the x-y plane and the unstable magnetic field component is thus 
aligned with the z-direction. Indeed only D^ z grows in our simulation, while the other 
magnetic field components remain at noise levels (not shown). 

Figure [6] shows the electron thermal anisotropy A(t) and compares the spatial 
distribution of B z at the time t = 150, which is just before the exponential growth 
phase of D\, z in Fig. E(a), and at t = 400, when the exponential growth phase of Db z 
ends. The magnetic field is initially weak and the structures are small. The saturation 
of the BTI between the bulk plasma and the electrons has resulted in a rapid growth of 
A, which has reached A ~ 1 at t — 150. This value corresponds to an electron thermal 
energy along x that is twice that along y. The anisotropy reaches its peak at t « 250 
and decreases thereafter to a value A »s 0.4 at t = 400, when a large coherent magnetic 
field structure has emerged in Fig. [6^b). The size of this patch is of the order of an 
electron skin depth in both spatial directions and the magnetic amplitude gives a ratio 
between the electron cyclotron frequency and the plasma frequency of 6 x 10 -4 , which 
is comparable to the ratio found in the interstellar medium [32J. The field patch is not 
thermal noise. Thermal noise [S] is linked to charge and current density fluctuations due 
to the finite number of computational particles and its coherence length is comparable 
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to the grid cell size. The structure centred at x = 0.7 and y = 1.7 in|^b) is, however, 
not the plane wave with a wave vector parallel to the y-axis, which we would expect 
from a TAWI [21]. A magnetic structure is growing while A is large, but we can not 
attribute it to an elementary instability such as a TAWI. 

3.2. Saturation of the instability between electrons and the fast ion beam 

The comparison of Figs. [1(a) and 0(b) demonstrates that the mean electric field (E x ) 
in Fig. [2(b) has accelerated the bulk of the electrons from v x — 1 to v x ~ 1.7. The 
velocity gap Av x between the value v x , where the electron phase space density reaches 
its maximum, and the mean speed Vj, = 6 of the fast beam has thus decreased from 
5 to 4.3. The change of Av x should result in a larger wave number of the unstable 
electrostatic waves than the now normalized k u2 = l/(vb — v e ) ~ 13 estimated in section 
2. This relative speed Av x expressed in units of the electron thermal speeds is also well 
below the value used in Ref . [H] , in particular if we take into account that the electrons 
have been heated up by the saturation of the first BTI. 

Figures E(a,b) reveal that this second BTI has already commenced to grow at 
t = 400. Movie 1 shows a wave perturbation after t w 300 of the electron distribution 
that propagates at a high speed to increasing x, which corresponds to the wave 
driven by the BTI between the electrons and the fast beam. Strong modulations 
of the electron distribution with 1.5 < v x < 3 are present. Their wave number is 
k a = 147r/3 ~ 14.6, which is indeed larger than k u2 . However, k a may not correspond 
to the true fastest growing wave because the periodic boundary conditions enforce a 
discrete wave spectrum. The electron velocity distribution in Fig. [7(c) reveals that the 
dilute electron component forms an almost circular distribution centred at v x ~ 0.5 and 
v y m 0. The contour line corresponding to 10 -1 ' 18 at v x ~ 1 is almost aligned with the 
Vy-axis. It resembles the crescent in Ref. [18], but it is much less pronounced. 

Figure [8] shows the electric field distribution at t — 500, when D ex and D ey reach 
their peak values in Fig. [2(a). The structures in E x and E y show a spatial correlation due 
to the obliquity of the wave front. The characteristic tilt angle between the wave vector 
and the x-axis is less than the one in Fig. [2] and thus \E y \ <C \E X \. Internal filamentary 
structures are visible in the intervals with \E X \ ^> 0. The filamentary structures are 
relatively stronger in E y compared to E x , which may explain why D ex does no longer 
grow in unison with D ey in Fig. [2(a) at this time. The electric field amplitude in Fig. 
[8(a) exceeds that in Fig. [3(a) by a factor of 4 and the wave length is larger by a factor 
5. The electrostatic potential, which is driven by the BTI between the electrons and 
the fast ion beam, thus exceeds that of the faster growing BTI by a factor of ~ 20 and 
should give trapped electron islands larger than those in Fig. [H 

Some electrons are trapped at the time t = 500 by this faster wave, which is the 
time when D ex saturates in Fig. [2(a). The supplementary movie 2, which shows the 
spatially integrated phase space distribution as a function of v x ,v y and corresponds to 
Fig. [9(c), demonstrates that the number of trapped electrons rapidly increases after this 
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Figure 7. (Color online) Electron phase space distributions at t = 400: Panel (a) and 
(b) show phase space projections onto the x, v x plane. The distribution is integrated 
over < y < 0.05 in (a) and over 1.55 < y < 1.6 in (b). The electron distribution has 
been integrated over all x,y,v z for the projection onto the v Xl v y plane in panel (c). 
All distributions are normalized to their respective maxima at t = 150 and the colour 
scale is 10-logarithmic. The black curves in (c) are the contour lines -1.2, -1.19, -1.18, 
-1.17, -1.16. 
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Figure 8. (Color online) The electric field distributions at t = 500: E x is shown in 
panel (a) and E y in panel (b). 



time. The black lines correspond to the speeds of both ion beams. Figures |9]Ja,b) show 
that the electron velocity distribution is not sinusoidal, which is a sign of a nonlinear 
wave, and that some electrons have speeds of v x > 5. These electrons are found close 
to the cusps of the electron distribution with v x ~ 5, which correspond to the unstable 
equilibrium point of a periodic electrostatic potential that moves with a speed ~ v^. The 
fastest electrons have just started their periodic motion in the electrostatic potential of 
the wave fronts with a wave vector, which is parallel to the x-axis. 

The electron velocity distribution in Fig. Ws) reveals that a small fraction of 
electrons have been accelerated up to v x ~ 10. These electrons have reached the stable 
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Figure 9. (Color online) Electron phase space distributions at i = 500: Panel (a) and 
(b) show phase space projections onto the x, v x plane. The distribution is integrated 
over < y < 0.05 in (a) and over 1.55 < y < 1.6 in (b). The electron distribution has 
been integrated over all x,y,v z for the projection onto the v x ,v y plane in panel (c). 
All distributions are normalized to their respective maxima at t = 150 and the colour 
scale is 10-logarithmic. 



equilibrium point of the electrostatic potential and their kinetic energy reaches its peak 
value. The electron distribution extends here and at later times (Movie 2) over a velocity 
interval that does not exceed 20v e or 0.3c and the electrons thus remain nonrelativistic. 
The velocity distribution reveals non-thermal electrons at v x ~ 1 and at \v y \ ~ 4. These 
electrons must have gained speed by their interaction with the oblique modes. They 
can interact more easily with electrons than the beam aligned mode, because their 
phase speed is lower [18] . The interaction of electrons with the oblique modes results 
in a velocity distribution, which resembles a crescent. This saturation mechanism is 
dominant, if the thermal spread of the electrons is small compared to the beam speed. 
The beam speed was 10 times the electron thermal speed in Ref. [18]. If the electron 
thermal spread is comparable to the beam speed, like for the instability between the bulk 
ions and the electrons, then the beam aligned modes can easily trap the electrons in the 
high-energy tail of the electron distribution and no crescent develops. The larger gap 
between the speed v x ~ 3.5 of the fastest electrons in Fig. 0(c) and the ion beam speed 
v b = 6 apparently puts the electron interaction with oblique modes and with beam- 
aligned modes on an equal footing. Movie 2 clearly shows that electrons are accelerated 
along v x at large v x and along v y for v x ~ after t w 450. 

The interaction of the electrons with the electrostatic potentials of the long waves 
scatters them. Scattering the electrons results in dissipation, which affects the current 
and charge density distribution of the plasma. Figure [10] shows the distributions of E x 
and B z at Tm = 1200. Strong electrostatic waves are still present. The magnetic field 
reveals strong noise fluctuations on a Debye-length scale, which we attribute to thermal 
noise. The typical amplitudes of B z are large; the fluctuation amplitude increases with 
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Figure 10. (Color online) The electric field distributions at t = 1200: E x is shown in 
panel (a) and B z in panel (b). 




Figure 11. (Color online) Particle distribution functions at t = 1200: Panel (a) 
shows the electron velocity distribution. Panel (b) and (c) show the ion distributions 
integrated over the interval < y < 0.05. All distributions are normalized to their 
peak value at t = 1200 and the colour scale is 10-logarithmic. 



the electron temperature. We conclude from Fig. ITOT b) that the TAWI, even if it has 
been responsible for the magnetic field growth in Fig. [6j is not an efficient source of 
coherent magnetic fields for the selected parameters. 

The reason for the discrepancy between the lifetimes of the electrostatic and 
magnetic structures becomes evident from the plasma distribution functions shown in 
Fig. HU The electron velocity distribution in Fig. HIT a) has its peak density at v x ~ 
and v y ~ 0. The electrons accumulate at speeds, which correspond to the phase speed of 
the waves driven by the BTI between bulk ions and electrons. An approximately circular 
electron velocity distribution extends up to ^Jv% + w 3. The electron scattering is 
isotropic in the rest frame of the slow waves. A second circular structure, albeit with 
a lower number density, is observed at v x ~ 6. These electrons accumulate at the 
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phase velocity of the waves driven by the BTI between electrons and the fast ion beam. 
The electron scattering is apparently isotropic in the rest frame of the fast waves. The 
electrons are thus scattered by two systems of electrostatic waves that have a different 
wave length and that move at different speeds. 

The scattering should result in the rapid thermalisation of the electrons and in 
a dissipation of large-scale electron currents in the simulation plane. This dissipation 
explains why the B z field in Fig. [TOT b) is no longer coherent over large spatial scales. 
The electrostatic fields can also be supported by modulations of the ion velocity and 
density and they are thus robust against the electron scattering. The beam ions in 
Fig. rrTTb) do not show any non-thermal features and the density modulation is in the 
linear regime. Their high temperature implies that the beam ions can support the large 
electrostatic fields in Fig. HUT a) without developing nonlinear velocity oscillations. The 
bulk ions in Fig. UTT c) show strong velocity oscillations. Their low thermal pressure 
implies that they can not sustain large electrostatic fields without showing non-thermal 
signatures, like accelerated ions. 

4. Discussion 

We have examined the growth and saturation of electrostatic Buneman-type instabilities 
(BTI's) with the help of a particle-in-cell simulation. Our system of two counter- 
streaming ion beams is representative for the foreshock region of supernova remnant 
shocks. The bulk ions and the electrons are the upstream plasma, while the ion beam 
corresponds to the shock-reflected ions. The system has initially been charge and 
current neutral and the plasma has been field-free, which is idealized. Electrostatic 
unmagnetized shocks have been observed experimentally and numerically |45l [46] , but 
their speed is typically well below the electron thermal speed. Their shock-reflected ion 
beam can thus not be faster than the electron thermal speed. We assume here that the 
SNR shock carries a magnetic field perpendicular to the shock normal, which reflects the 
incoming upstream ions even at high shock speeds, that is weak enough to be negligible 
with regard to the growth of the BTI. The latter implies that, in the electron rest frame, 
the electron cyclotron frequency should be well below the frequency of the waves driven 
by the BTI. A stronger magnetic field yields additional unstable wave branches [16J and 
isotropizes the electrons more rapidly orthogonally to the magnetic field. 

The simulation has shown that two wave modes are growing. One instability branch 
corresponds to the BTI between the bulk ions and the electrons. It is the faster growing 
one. The second branch develops between the fast ion beam and the electrons. The large 
phase speed of its waves can accelerate electrons to larger speeds and this instability 
is thus more powerful. The electron speeds do, however, remain nonrelativistic as in 
previous ID simulation studies [13]. The BTI in unmagnetized plasma can thus not 
account for electron injection into the diffusive acceleration and magnetic fields would be 
necessary for this purpose [47J. Both wave branches support a continuos wave spectrum 
ranging from the beam-aligned modes with a high phase speed to the slower oblique 
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modes [TS]. The low relative speed between the electrons and the bulk ions implied 
that the instability saturated by the conventional trapping of electrons by beam-aligned 
modes [34|. The growth of magnetic fields could be observed after this saturation. 
Their energy density was, however, too low to identify unambigously the instability 
responsible for their growth. A likely candidate is the TAWI [31]. The peak ratio 
between the electron cyclotron frequency and the plasma frequency remained below 
10 -3 , which is comparable to the equivalent in the ISM. 

The electrostatic instability between the fast ion beam and the heated electrons 
developed at a later time. It saturated by the trapping of electrons by beam-aligned 
modes and by obliquely propagating slower modes. The phase space distribution 
revealed the simultaneous development of phase space holes and of the crescent 
distribution, which was reported first in Ref. [18]. The second mechanism becomes 
dominant, if the beam speed exceeds by far the electron thermal spread, which was not 
the case here. The formation of such a distribution might also be aggravated by the 
much hotter ions we have used here. No growth of a strong coherent magnetic field 
could be observed, in spite of a stable thermal anisotropy of the electron distribution. 

It is unclear, why the TAWI is inefficient here. A possible reason is that the 
magnetic interaction between electrons, which gives rise to the TAWI, competes with 
the electron interaction with the electrostatic wave fields. The latter has not been 
taken into account in numerical studies and may inhibit the formation of stable current 
filaments that are needed to sustain strong magnetic fields. We conclude that, at least 
for our simulation parameters, the TAWI driven by the BTI can not magnetize SNR 
shocks due to its short life time and low peak amplitude of the magnetic fields. 
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